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Abstract 

In contrast with most Monte Carlo integration algorithms, which are used to estimate ratios, the 
replica gas identities recently introduced by Adib enable the estimation of absolute integrals and 
partition functions using multiple copies of a system and normalized transition functions. Here, an 
optimized form is presented. After generalizing a replica gas identity with an arbitrary weighting 
function, we obtain a functional form that has the minimal asymptotic variance for samples from 
two replicas and is provably good for a larger number. This equation is demonstrated to improve 
the convergence of partition function estimates in a 2D Ising model. 
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The standard paradigm in Monte Carlo integration is to estimate a ratio. In the prototyp- 
ical scenario, a simulation is used to obtain the fraction of randomly sampled points which 
lie within an encompassing reference region. The integrand of interest is then estimated by 
multiplying this fraction with the known reference area. 

More sophisticated algorithms retain the feature that they estimate ratios of normalizing 
constants [l| or partition functions |2|. Furthermore, the efficiency of all these integration 
techniques depend on the degree of overlap between the target integrands. Hence, much 
effort has been devoted to obtaining intermediate functions which bridge the integrands of 
interest. 

Recently, Adib introduced a class of theorems, which he dubbed the replica gas identities, 
that can be used to estimate absolute integrals [3;|. These identities are based on the simple 
physical idea of measuring the volume of a container by filling it with a specified density 
of ideal gas, and then counting the number of particles. Similarly, an integrand can be 
estimated by creating multiple non-interacting copies (or replicas) of a system that fill a 
region with a specified density, and then counting the number of replicas. Adib also described 
a complementary and more abstract identity relevant to simulations with a fixed number of 
replicas. A variant of this latter identity, which is particularly applicable to existing parallel 
tempering (a.k.a. replica exchange) simulations, will be the focus of the present work. 

In this work, we shall consider integrals of the form, 



where Q is the support of the integral and tt(x) is a positive-definite function of a 
d— dimensional vector x (see Ref. for a generalization to non-positive definite integrands). 
The function tc(x) can be regarded as an unnormalized probability density. For many equili- 
brated physical systems with the energy function E(x), %(x) = e~ E ^ describes the relative 
probability of observing x, and Z is the partition function. Discrete integrals can be treated 
analogously. 

Suppose that we have multiple replicas of a system on the same support, each inde- 
pendently sampling from their own respective distributions. The replica gas identities cou- 
ple these copies together by means of transition functions, T(x'\x), which are normalized 
conditional probability distributions that one can sample from and evaluate. Appropriate 
transition functions include candidate-generating functions routinely used in Markov Chain 
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Monte Carlo simulations, such as a Gaussian centered about x. 

For simplicity, we shall first treat the case where the number of replicas is fixed at two. 
One replica samples from tt(x) and the other samples from an auxiliary distribution tt(x) 
(which has the auxiliary integral Z = J n dx7r(x)). The density tt(x) is arbitrary and may 
be ir(x) itself, but usually corresponds to ir(x) at a different temperature. In this case, a 
generalization of the second replica gas identity [3] is, 

j n dx [tc(x)/Z] J n dx' T(x'\x) ■ a(x', x) ir(x') 



f n dx [tt(x)/Z] f n dx' [ir(x')/Z] ■ a(x',x) T(x'\x) 
(air(x'))~ T 



(2) 



(aT{x'\x))~^ 

where the expectation (0(x, x'))j g of an observable 0(x, x') means that x and x' are sampled 
from the distributions / and g, respectively. In the numerator of the equation above, x is 
sampled from 7r(x), x' is sampled using transition function, and tc(x') is evaluated. In the 
denominator, x is also sampled from tt(x), but x' is drawn from the replica with density 
7r(x'), and T[x'\x) is evaluated. The identity also includes an arbitrary weighting function 
a = a(x',x); any function is valid as long as the denominator is nonzero. The choice a = 1 
reduces the expression to Adib's original identity. 

The replica gas identity can be made into an estimator by replacing the expectations 
with sample means. While many choices of a will yield statistically consistent expressions, 
some choices will lead to more efficient estimators than others. Here, we will optimize a 
by minimizing the asymptotic variance in estimates of C = InZ, which we shall denote by 
(. Our procedure will closely follow that of Bennett [4|, who minimized the error in an 
estimator for the ratio of partition functions, and yield an analogous estimator. 

In the large sample limit, ( will be Gaussian about the true value of (. Using error 
propagation based on a first-order Taylor series expansion, the asymptotic variance is, 

aT] = («M^)%,t , (a 2 T{x'\xY)^ 1 1 
N T (a7r(x'))l T N^aTix'lx)) 2 ^ N T N v 
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(aT(x'\x))l n N T N n ' 

where N T and N w are the number of independent samples used to estimate (0(x,x'))~ T 
and (0(x, x'))~ 7r , respectively. The second line is obtained by writing both expectations in 
terms of (0(x, x'))~ n and combining terms. 



As multiplying a by a constant does not change the value of <J 2 [Q, we can use a constraint 
where {aT(x'\x))- n (and hence the denominator) is constant. The variance, Eq. [3], is then 
minimized using Lagrange multipliers, leading to the optimal a, 



(4) 



Notably, a* includes the sought quantity (. The estimator obtained by substituting this 

function into Eq. [2] can either be solved by self-consistent iteration, as originally described 
I I 

by Bennett [4|], or by finding the zero of an implicit function obtained by rearrangement, 
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which is similar to an expression described by Shirts et. al. [5]. The variance in ( in Eq. [5] 
may be estimated by substituting the appropriate a into Eq. [3j 
By substituting a* into Eq. [3], multiplying by a factor of unity, 
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and separating the expression into two expectations, we obtain the following formula for the 
variance, 
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which is analogous to an extant expression for the error in the Bennett's method 

When samples are drawn from K replicas with integrands TTk(x) and partition functions 
Zfc for k = 1, K, the following generalized replica gas identity is applicable: 



(8) 



Here, 7r t (:r) and Z t are the targeted integrands and integrals, respectively. The sums run 
over all replicas which are not t. As in Eq. [2J the choice a^ — l yields an identity originally 
described by Adib {j}]. 

We now consider the optimization of a in the case where K > 2. The asymptotic variance 
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of Ct is, 
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Unlike in the case where K = 2, the second term in the above variance expression does not 
reduce to a constant. We can show, however, using a procedure similar to Guibas and Veach 
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7J, that it is bounded by a constant: 
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Thus, if we find a set of functions for en*, that minimize the first term in cr[Ct], no other 
choices for at can decrease the variance by more than max ( — h i^— ) . 

Minimization of the first term in Eq. [9] by Lagrange multipliers leads to an optimized set 
of a k analogous to Eq. HI As with K = 2, substitution of this a k leads to an estimator for 
the integral which can either be solved self-consistently or by finding the zero of an implicit 
function, 
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The corresponding variance may be estimated by substituting the appropriate ai^ into Eq. 
9. 

The original and optimized forms of the replica gas identity were demonstrated by esti- 
mating the partition function of a two-dimensional 32 x 32 spin Ising model, using parame- 
ters similar to those previously described fl. As in a typical Ising model simulation without 



magnetization, the energy of a spin configuration x was given by E(x) = — J2<k i> x k x h a 
sum over all pairs of neighbors 8]. Periodic boundary conditions were used. Simulations 
were performed for 26 replicas at the temperatures = 0.5, 0.6, 3.0. Configurations 
were equilibrated using 250 iterations in which a cluster move, Woolf's generalization {9] 
of the Swendsen-Wang algorithm |lo[, was followed by attempted replica exchanges be- 
tween adjacent temperature pairs (alternating between #2), (/?3, At), (Pk-i, Pk)} an d 
{(/?2j ^3), (/?4, ^5), {Pk-2, Pk-i)})- Equilibration was followed by 10 5 production steps us- 
ing Monte Carlo trial moves generated by a transition function T(x'\x) that flips each spin 
with probability Pfu p = 1/N, where N is the total number of spins. The transition function 
was the same for every replica. After every 100 steps of production, a replica exchange was 
attempted between a random replica and its higher-temperature neighbor, energies of each 
configuration were stored, T(x'\x) was used to sample configurations from (0(x, x')) n T 
(these are distinct from the Markov chain), and transition probabilities between the dif- 
ferent configurations were calculated using T{x'\x) = (p/ii P )" x '~ x "(l — Pfn p ) N ~^ x ~ x ^, where 
\\x' — x\\ is the number of spin mismatches 3j. 

For this system, we find that both estimators have similar and relatively accurate ef- 
ficiencies at lower temperatures (Fig. [I]). Closer to and above the critical temperature 
{ksT = 2.269), however, the optimized estimator (Eq. fT2!) has less variance and bias than 
the original estimator (Eq. [8]) with unit weight. Nonetheless, the high-temperature estimates 
are substantially worse than those at low temperatures. 

Ultimately, the quality of the optimized expression is subject to the same fundamen- 
tal limitations as Adib described with the original estimator 3[]. For the convergence of 
averages to not be dominated by rare events, (a) most configurations generated by the tran- 
sition function should fall in typical regions of phase space for the target distribution and 
(b) most pairs of configurations from different replicas should have substantial density in 
T(x'\x). These requirements hold true in the low-temperature regime, where the spins are 
relatively ordered and a local transition function suffices, but are not met at higher temper- 
atures, where configurations have disordered spins. As Adib mentioned, convergence may 
be improved by fine-tuning the simulation procedure or choosing a more suitable transition 
function. 

Overcoming this sampling limitation is not the intent of the current work. Rather, the 
purpose is to improve the analysis of a given data set. While the resulting expression (Eq. 



6 



5000 



N 



4000 
3000 
2000 
1000 




0.5 



2.5 



p-(=k B T) 



FIG. 1: Comparison of logarithmic partition function estimates: the exact Kaufman formula 
(dashed line), the original replica gas identity, Eq. [8] with a = 1 (upward triangles), and the 
optimized form, Eq. [12] (downward triangles). Markers and error bars indicate the mean and 
standard deviation of 1000 independent simulations. Inset: magnification around the critical tem- 
perature. 



[T2|) is more complex, the increased difficulty of implementation is no greater than using the 
Bennett Acceptance Ratio [4] relative to unidirectional estimates of free energy differences, 
and does not require sampling from more ensembles. Expected statistical efficiency gains 
should make this investment worthwhile. 
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